# Generate a buffer with 2km on either side of a to-code / no-code border (i.e., either SRA or LRA-H on one side, LRA not H on the other)
pacman::p_load(tidyverse, mapview, sf, tigris, qs, geos)
sf_use_s2(FALSE)

source("code/globals.R")

# Load SRA areas
responsibility_areas <- st_read(file.path(INPUT_PUBLIC, "SRAandFHSZ/State_Responsibility_Area-shp/Responsibility_Areas.shp"))
responsibility_areas <- responsibility_areas %>% st_transform(3310)

SRAonly <- responsibility_areas %>% filter(SRA == "SRA") %>% st_union() %>% st_sf(SRA = "SRA")
LRAonly <- responsibility_areas %>% filter(SRA == "LRA") %>% st_union() %>% st_sf(SRA = "LRA")

# Load LRA Hazard classifications --- Ideally these should be the same year to avoid weird issues
hazard_class <- st_read(file.path(INPUT_PUBLIC, "SRAandFHSZ/LRA_FHSZ/California_Fire_Hazard_Severity_Zones__FHSZ_.shp")) %>%
  st_transform(3310)

LRAHonly <- hazard_class %>%
  filter(complete.cases(SRA, HAZ_CLASS)) %>%
  transmute(LRA_VHFHSZ = TRUE) %>%
  st_union() %>% 
  st_sf(LRA_VHFHSZ = TRUE)

# Construct unioned tocode and nocode areas
tocode <- SRAonly %>% st_union(LRAHonly) # to code is SRA plus LRA-VHFHSZ
nocode <- LRAonly %>% st_difference(LRAHonly) # no code is LRA minus LRA-VHFHSZ

# Shrink and expand buffer by 100m to get rid of slivers
nocode <- nocode %>% st_buffer(-100) %>% st_buffer(100)

compute_border <- function(border_size = 1000) {
  print(sprintf("Computing border for %s meter buffer", border_size))
  # Create a polygon that is the border of SRA/LRA regions
  # border_size is in meters

  # Buffer by half the desired border size
  buffer_tocode <- tocode %>% st_buffer(border_size / 2)
  buffer_nocode <- nocode %>% st_buffer(border_size / 2)

  # Compute intersection between SRA/LRA buffers
  border <- st_intersection(buffer_tocode, buffer_nocode)

  # Return the union
  border %>% st_union()
}

border_250m <- compute_border(250)
border_500m <- compute_border(500)
border_1000m <- compute_border(1000)
border_2000m <- compute_border(2000)

# Create an SF object with borders
borders <- st_sf(
  border_size = c(250, 500, 1000, 2000),
  geometry = c(border_250m, border_500m, border_1000m, border_2000m)
)

new_borders <- borders
qsave(borders, file.path(WORKING, "sra-lra-borders.qs"))

